#include "ann_math.h"
#include "tnt/tnt.h"
#include "tnt/jama_svd.h"

#include <iostream>

namespace ann
{

namespace math
{
    void transpose( TNT::Array2D<double>& A, TNT::Array2D<double>& At )
    {
        UINT m = A.dim1();
        UINT n = A.dim2();

        for( UINT i = 0; i < m; ++i )
        {
            for( UINT j = 0; j < n; ++j )
            {
                At[j][i] = A[i][j];
            }
        }
    }

    const UINT SolveGwt( double* pG, const UINT pm, const UINT pn, double* pt, double* pw )
    {
        TNT::Array2D<double> G(pm, pn, pG);
        TNT::Array2D<double> Gp(pn, pm, 0.0);
        TNT::Array2D<double> U(pm, pn);
        TNT::Array2D<double> Ut(pn,pm);
        TNT::Array2D<double> S(pn, pn);
        TNT::Array2D<double> V(pn, pn);
        TNT::Array1D<double> t(pm, pt);
        TNT::Array1D<double> w(pn, 0.0);

        JAMA::SVD<double> svd(G);
        svd.getU(U);
        svd.getS(S);
        svd.getV(V);

        ann::math::transpose(U, Ut);
        for (UINT i = 0; i < pn; ++i)
        {
            S[i][i] = 1.0 / S[i][i];
        }
        Gp = TNT::matmult(TNT::matmult(V, S), Ut);

        std::fill_n(pw, pn, 0.0);
        for( UINT i = 0; i < pn; ++i )
        {
            for( UINT j = 0; j < pm; ++j )
            {
                pw[i] += Gp[i][j]*pt[j];
            }
        }

        return 0;
    }
}

}
